function [ratio_trimmed_mean,mean_ratio,trimmed_mean_ratio] = script_ratio_rand_exp_fun_average(theta_NPD,m_regr,data_wide,J,ntheta,combos_all_unique)

N_temp = 500;

% combos_all_unique = combos_all;
combos_all_unique = double(combos_all_unique);

ratio_trimmed_mean_temp = NaN(2,1);
mean_ratio_temp = NaN(2,1);
trimmed_mean_ratio_temp = NaN(2,1);

for j=1:2
    
    ratio_temp = zeros(N_temp,1);
    d2s_dzdz = zeros(N_temp,1);
    d2s_dzdx = zeros(N_temp,1);
    der_num = zeros(N_temp,ntheta);
    der_den = zeros(N_temp,ntheta);
    
    parfor i=1:N_temp
        
        [~,temp] = max(data_wide.z(i,:));
        
        temp_diff = setdiff(1:J,temp);
        
        if j==2
            
            temp_diff = temp_diff([2 1]);
            
        end
        
        ind_good2 = temp_diff(1);
        
        arg_x_1 = normcdf((data_wide.x(i,temp)-mean(data_wide.x(:)))/std(data_wide.x(:)));
        arg_x_other = normcdf((data_wide.x(i,temp_diff)-mean(data_wide.x(:)))/std(data_wide.x(:)));
        
        arg_z_1 = normcdf((data_wide.z(i,temp)-mean(data_wide.z(:)))/std(data_wide.z(:)));
        arg_z_other = normcdf((data_wide.z(i,temp_diff)-mean(data_wide.z(:)))/std(data_wide.z(:)));
        
        argument = [arg_x_1 arg_x_other arg_z_1 arg_z_other];
        
        deriv_z_1 = normpdf((data_wide.z(i,temp)-mean(data_wide.z(:)))/std(data_wide.z(:)))*1/std(data_wide.z(:));
        deriv_z_2 = normpdf((data_wide.z(i,ind_good2)-mean(data_wide.z(:)))/std(data_wide.z(:)))*1/std(data_wide.z(:));
        deriv_x_2 = normpdf((data_wide.x(i,ind_good2)-mean(data_wide.x(:)))/std(data_wide.x(:)))*1/std(data_wide.x(:));
        
        [temp_num,temp_der_num] = fun_d2sigma_new2(theta_NPD,argument,m_regr,combos_all_unique,J,J+1,J+2);
        d2s_dzdz(i) = temp_num*deriv_z_1*deriv_z_2;
        der_num(i,:) = temp_der_num*deriv_z_1*deriv_z_2;
        
        [temp_den,temp_der_den] = fun_d2sigma_new2(theta_NPD,argument,m_regr,combos_all_unique,J,J+1,2);
        d2s_dzdx(i) = temp_den*deriv_z_1*deriv_x_2;
        der_den(i,:) = temp_der_den*deriv_z_1*deriv_x_2;
        
        ratio_temp(i) = d2s_dzdz(i)/d2s_dzdx(i);
        %i
    end
    
   [num_temp_ord,ind_num] = sort(d2s_dzdz);
    [den_temp_ord,ind_den] = sort(d2s_dzdx);
    index = floor(N_temp*0.25):floor(N_temp*0.75);
    ratio_trimmed_mean_temp(j) = mean(num_temp_ord(index))/mean(den_temp_ord(index));
    
  
    mean_ratio_temp(j) = mean(ratio_temp);

    [temp_ord,ind] = sort(ratio_temp);
    index = floor(N_temp*0.25):floor(N_temp*0.75);
    trimmed_mean_ratio_temp(j) = mean(temp_ord(index));

end

ratio_trimmed_mean = mean(ratio_trimmed_mean_temp);

mean_ratio = mean(mean_ratio_temp);
trimmed_mean_ratio = mean(trimmed_mean_ratio_temp);

end
